Análisis Exploratorio de Suelos (Honoria - Tournavista)

Publicación Internacional:

Certificado de Aceptacion:

Study to Determine Levels of Cadmiumin Cocoa Crops Applied to Inland Areas of Peru

Paper:

Acreditation

Link de Descarga

El link es el siguiente: Agronomy MDPI.

Video de Presentación

El link es el siguiente: ECITEC UNI 2020

Análisis de la Información Recolectada:

Revisión del Modelo de Elevación Digital:

DEM <- raster("Data_Suelos/raster/DEM_DRONE/Zona_Oeste.tif")
# str(DEM)
# DEM
# slope <- terrain(DEM, opt ="slope")
# aspect <- terrain(DEM, opt="aspect")
# DEM.hill <- hillShade(slope, aspect, angle = 40, direction = 270)
DEM <- crop(DEM, extent(518600, 519200, 9031800, 9032400))
#plot(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
#     zlim=c(0, 250))
image(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
      zlim=c(0, 250))

Revisión del DEM con factores:

Factores_cacao <- st_read(dsn = "Data_Suelos/shp/Cacao_Factors/Age_Planta.shp")

## Reading layer `Age_Planta' from data source 
##   `D:\CursoDAG\DataAnalisisGeociencias\Modulos\ModuloII\Data_Suelos\shp\Cacao_Factors\Age_Planta.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 518541 ymin: 9031841 xmax: 519584 ymax: 9032609
## Projected CRS: WGS 84 / UTM zone 18S

mapview(Factores_cacao)

mapview(DEM, legend =TRUE)+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
         lwd = 2, color = "blue")

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 573840800 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  573840800 '

Generando Shapefiles (Opcional):

Factores_cacao2 <- as.data.frame(Factores_cacao)
#Factores_cacao2$Uso_Tierra <- edit(Factores_cacao2$Uso_Tierra)

Factores_cacao2$Var_cocoa <- c("Creole-Aromatic","CCN51","","CCN51","Creole-Aromatic",
                               "CCN51")
Factores_cacao3 <- st_as_sf(Factores_cacao2)

mapview(Factores_cacao3, zcol="Var_cocoa")

Honoria_elementos <- readxl::read_xlsx(path = "Data_Suelos/csv_xlsx_txt/BD_Cd.xlsx")
str(Honoria_elementos)

## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Xeast         : num [1:24] 518642 518644 518690 518701 518706 ...
##  $ Ynorth        : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
##  $ Cdsoil_mgkg   : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
##  $ Pbsoil_mgkg   : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
##  $ Znsoil_mgkg   : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
##  $ pHapprox      : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
##  $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...

Honoria_elementos <- Honoria_elementos %>% rename(Este= Xeast, Norte=Ynorth, 
                                                  Cd = Cdsoil_mgkg,
                                                  Pb = Pbsoil_mgkg,
                                                  Zn = Znsoil_mgkg,
                                                  pH = pHapprox,
                                                  Ce = Ceapprox_uS_cm)
Honoria_elementos2 <-Honoria_elementos[ ,c("Norte","Este")]
Honoria_elementos2 <- Honoria_elementos2[ ,order(c(names(Honoria_elementos2)))]
sputm  <- SpatialPoints(Honoria_elementos2, proj4string=CRS("+proj=utm +zone=18 +south +datum=WGS84")) 
spgeo  <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
spgeo  <- as.data.frame(spgeo)
colnames(spgeo) <- c("lng","lat")
Honoria_elementos2 <- cbind(Honoria_elementos,spgeo)

Honoria_elementos3 <- st_as_sf(Honoria_elementos2, coords = c("lng", "lat"),
                               remove = FALSE, crs = 4326, agr = "constant")

mapview(Honoria_elementos3, zcol = "Cd") +mapview(Factores_cacao3, zcol="Var_cocoa")+
  mapview(Factores_cacao3, zcol = "Uso_Tierra")+
  mapview(Factores_cacao3, zcol = "Age")

Factores_cacao4 <- Factores_cacao3 %>%
  st_transform(crs = 4326 )
  
Honoria_elementos3 %>%
  st_intersection(Factores_cacao4) %>%
  mapview(zcol = "Uso_Tierra", burst = TRUE, cex = "Cd")+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")

Honoria_elementos3 %>%
  st_intersection(Factores_cacao4) %>%
  mapview(zcol = "Age", burst = TRUE, cex = "Cd" )+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")

Honoria_elementos3 %>%
  st_intersection(Factores_cacao4) %>%
  mapview(zcol = "Var_cocoa", burst = TRUE, cex = "Cd")+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")

Estandar de Calidad de Cd en Suelos:

\[Cd<=1.4\]

mapview(Honoria_elementos3, zcol = "Cd", at = c(0,1.4,2.5), legend = TRUE,
         col.regions = c("green", "red"))+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")

Data Table:

Factores_cacao_n <- Factores_cacao 
Honoria_elementos_n <- Honoria_elementos3 %>% st_transform(crs = st_crs(Factores_cacao))
# sf_use_s2(TRUE) #cambiar a spherical geometry
intersection <- st_intersection(Honoria_elementos_n, Factores_cacao_n)
intersection

Honoria_statical <- intersection %>% st_set_geometry(NULL)

Statical Summary:

min_max_mean_sd <- list(
  min = ~min(.x, na.rm = TRUE),
  max = ~max(.x, na.rm = TRUE),
  mean = ~mean(.x, na.rm = TRUE),
  sd = ~sd(.x, na.rm = TRUE)
)

Honoria_statical %>% 
  summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))

Honoria_statical %>% group_by(Age) %>% 
  summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))

Ver vignette("colwise") para más detalles.

Boxplot Analysis:

Cd total:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  coord_flip()

Cd by Factor Age and Land Use:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Age)+
  coord_flip()

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Uso_Tierra)+
  coord_flip()

Cd by Prior Exploration Field:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Uso_Tierra)+
  coord_flip()

Average Cadmiunm Soil

Corrplot Bivariant

correlacion <- cor(Honoria_statical[ ,c("Cd", "Pb", "Zn", "pH", "Ce")]
                   , method = "pearson")
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria", 
         hclust.method = "median", addrect = 2)
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria")

correlacion2 <- cor(Honoria_statical[ Honoria_statical$Uso_Tierra=="Maiz 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
# correlacion2 <- cor(Honoria_statical %>% filter(Uso_Tierra=="Maiz 10") %>%
#                     select(Cd, Pb, Zn, pH, Ce), method = "pearson")
corrplot(correlacion2, method ="number")

corrplot(correlacion2, method ="circle", order = "AOE")

corrplot(correlacion2, method = 'shade', order = 'AOE', diag = FALSE) 

corrplot(correlacion2, method = 'square', order = 'AOE', diag = FALSE, type = "lower")

correlacion3 <- cor(Honoria_statical[Honoria_statical$Uso_Tierra=="Pasto 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
corrplot(correlacion3, method ="circle")

corrplot.mixed(correlacion3, order = 'AOE')

Para ver toda la funcionalidad revisar : An Introduction to corrplot Package.

Biplot Multivariant (Principal Component Analysis)

Preparacion de la informacion

library(FactoMineR); library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Honoria_pca_var <- Honoria_statical 
Honoria_pca <- Honoria_pca_var %>% dplyr::select(Cd, Pb, Zn, pH, Ce)

Usando la funcion prcomp()

res.pca <- prcomp(Honoria_pca, scale = TRUE)
res.pca <- PCA(Honoria_pca, graph=FALSE)
print(res.pca)

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 22 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

res.pca$call$centre

## [1]  1.664091  5.980000 59.403182  6.530000 27.547727

Eigenvalores

# La varianza explicada por el primer eigenvalor es 38.60% y el segundo 25.58% los cuales
# seran tomados para el analisis de los componentes principales representando 64.19% de la
# varianza total. Según (Kaiser 1961) un eigenvalue>1 indica un punto de corte 
# para contener el deja a voluntad del investigador realizar el corte en un valor adecuado
# para la visualizacion. 
eig.val<- get_eigenvalue(res.pca)
eig.val

##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1  1.9300862        38.601723                    38.60172
## Dim.2  1.2794547        25.589095                    64.19082
## Dim.3  0.9204006        18.408013                    82.59883
## Dim.4  0.6054406        12.108812                    94.70764
## Dim.5  0.2646179         5.292357                   100.00000

fviz_eig(res.pca, addlabels = TRUE, ylim=c(0,50))

Circulo de correlacion

var <-get_pca_var(res.pca)
var

## Principal Component Analysis Results for variables
##  ===================================================
##   Name       Description                                    
## 1 "$coord"   "Coordinates for the variables"                
## 2 "$cor"     "Correlations between variables and dimensions"
## 3 "$cos2"    "Cos2 for the variables"                       
## 4 "$contrib" "contributions of the variables"

#Circulo de Correlacion:
#La correlacion entre las variables de los componentes principales son usadas como
# coordenadas de la variable para los PC. Segun esto las observaciones son representadas
# por su proyeccion, pero las variables son representadas por sus correlaciones.
var$coord

##         Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
## Cd  0.5097511 -0.5694786  0.50836826 -0.35210702  0.18283969
## Pb  0.6597640 -0.4324834  0.04309302  0.61090780 -0.05103086
## Zn -0.2583970  0.5477468  0.75441302  0.24677884  0.05626362
## pH  0.7076697  0.5604771 -0.27468642  0.02387721  0.33022158
## Ce  0.8169370  0.3923593  0.12455405 -0.21629371 -0.34113264

fviz_pca_var(res.pca, col.var="black")

fviz_pca_var(res.pca, col.var = "cos2",
             gradient.cols=c("#00AFBB","#E7B800","#FC4E07"),
             repel=TRUE)

# El plot de correlacion nos indica que las variables pH y Ce estan correlacionadas
# positivamente (agrupadas), asi como el Cd, Pb. Las variables Cd-Zn y Pb-Zn tienen
# una correlación inversa, mientras que entre ph-Zn y Ce-Zn no existe una relacion
# significativa.

# Las variables pH, Ce, Pb y Cd estan bien representadas, adicionalmente
# se debe notar que el Zn disminuye su representatividad.

Contribucion de PCA

fviz_cos2(res.pca, choice="var",axes=1:2, top=10)

fviz_contrib(res.pca, choice="var",axis=1,top=10)

fviz_contrib(res.pca, choice="var",axis=2,top=10)

# La representatividad de la proyeccion para el Zn no es considerablemente
# significativa para las 
# Dimensiones 1 y 2, porque el angulo de proyeccion es muy alto por lo cual no 
# se proyecta de la manera deseada.

Correlacion y p-valores

res.desc <- dimdesc(res.pca, axes=c(1,2),proba=0.05)
res.desc$Dim.1

## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##    correlation      p.value
## Ce   0.8169370 3.498136e-06
## pH   0.7076697 2.296007e-04
## Pb   0.6597640 8.356745e-04
## Cd   0.5097511 1.537322e-02

res.desc$Dim.2

## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##    correlation     p.value
## pH   0.5604771 0.006664419
## Zn   0.5477468 0.008319053
## Pb  -0.4324834 0.044405255
## Cd  -0.5694786 0.005667244

Indicadores y Graficos PCA:

ind <- get_pca_ind(res.pca)
ind

## Principal Component Analysis Results for individuals
##  ===================================================
##   Name       Description                       
## 1 "$coord"   "Coordinates for the individuals" 
## 2 "$cos2"    "Cos2 for the individuals"        
## 3 "$contrib" "contributions of the individuals"

fviz_pca_ind(res.pca, col.ind="cos2",
             gradient.cols=c("#00AFBB","#E7B800","#FC4E07"),
             repel=TRUE)

fviz_contrib(res.pca, choice="ind", axes=1:2) #Contribucion puntual

fviz_pca_ind(res.pca,
             geom.ind = "point",
             col.ind = Honoria_pca_var$Age,
             palette = rainbow(6),
             addEllipses = TRUE,
             legend.title="Age")

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

fviz_pca_ind(res.pca,
             geom.ind = "point",
             col.ind = Honoria_pca_var$Uso_Tierra,
             palette = rainbow(4),
             addEllipses = TRUE,
             legend.title="Land Use")

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

q <- fviz_pca_biplot(res.pca, 
                   
                   geom.ind ="point",
                   fill.ind = Honoria_pca_var$Age, col.ind = "black",
                   pointshape = 21, pointsize = 2,arrowsize=1,
                   palette = rainbow(6),
                   addEllipses = TRUE,
                   ellipse.level= 0.95,
                   
                   
                   col.var = "black",
                   gradient.cols = "RdYlBu",
                   legend.title = "Time of \n Age",
                   repel=FALSE
)
q

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plotly::ggplotly(q)

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

Clusterizar la Data

# Compute PCA with ncp = 3
res.pca1 <- PCA(Honoria_pca, ncp = 6, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc <- HCPC(res.pca1, graph = FALSE)

Dendongrama

fviz_dend(res.hcpc, 
          cex = 0.7,                     # Label size
          palette = "jco",               # Color palette see ?ggpubr::ggpar
          rect = TRUE, rect_fill = TRUE, # Add rectangle around groups
          rect_border = "jco",           # Rectangle color
          labels_track_height = 0.8      # Augment the room for labels
          )

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

fviz_cluster(res.hcpc,
             repel = TRUE,            # Avoid label overlapping
             show.clust.cent = TRUE, # Show cluster centers
             palette = "jco",         # Color palette see ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map"
             )

Agregando informacion de la clusterizacion y visualizando

res.hcpc$desc.var$quanti

## $`1`
##       v.test Mean in category Overall mean sd in category Overall sd
## Ce -2.400916        14.236667    27.547727      9.9841602 27.8314601
## Cd -2.693181         1.470833     1.664091      0.2436342  0.3602229
## Pb -3.566253         5.405833     5.980000      0.5037933  0.8082135
##        p.value
## Ce 0.016354097
## Cd 0.007077392
## Pb 0.000362122
## 
## $`2`
##       v.test Mean in category Overall mean sd in category Overall sd
## Pb  3.700797         6.764444     5.980000      0.4474399  0.8082135
## Cd  2.685257         1.917778     1.664091      0.3517821  0.3602229
## Zn -2.243528        45.891111    59.403182     16.3646560 22.9641130
##        p.value
## Pb 0.000214923
## Cd 0.007247400
## Zn 0.024862798
## 
## $`3`
##      v.test Mean in category Overall mean sd in category Overall sd
## Ce 3.562238           126.69     27.54773              0 27.8314601
## pH 2.377132             8.00      6.53000              0  0.6183923
##         p.value
## Ce 0.0003677074
## pH 0.0174478606

Honoria_cluster <- res.hcpc$data.clust
View(Honoria_cluster)
library(tidyverse)
Honoria_cluster %>%
  filter(clust=="1")

Honoria_mapa_cluster <- cbind(Honoria_pca_var, Honoria_cluster[,c("clust")])
Honoria_mapa_cluster <- Honoria_mapa_cluster %>% dplyr::rename(clust = "Honoria_cluster[, c(\"clust\")]")

Honoria_mapa_cluster <- st_as_sf(Honoria_mapa_cluster, coords = c("lng", "lat"),
                               remove = FALSE, crs = 4326, agr = "constant")
mapview(Honoria_mapa_cluster, zcol="clust")

Cadmium Final Geostatistics Map:

library(gstat)

class(Honoria_elementos3)

## [1] "sf"         "data.frame"

summary(Honoria_elementos3)

##       Este            Norte               Cd              Pb       
##  Min.   :518642   Min.   :9031868   Min.   :1.100   Min.   :4.740  
##  1st Qu.:519033   1st Qu.:9031913   1st Qu.:1.350   1st Qu.:5.247  
##  Median :519183   Median :9032036   Median :1.575   Median :5.850  
##  Mean   :519108   Mean   :9032052   Mean   :1.634   Mean   :5.929  
##  3rd Qu.:519221   3rd Qu.:9032134   3rd Qu.:1.863   3rd Qu.:6.525  
##  Max.   :519556   Max.   :9032343   Max.   :2.550   Max.   :7.240  
##        Zn              pH              Ce              lng        
##  Min.   :11.80   Min.   :5.410   Min.   :  4.50   Min.   :-74.83  
##  1st Qu.:44.79   1st Qu.:6.160   1st Qu.:  9.00   1st Qu.:-74.83  
##  Median :58.70   Median :6.375   Median : 13.50   Median :-74.83  
##  Mean   :58.66   Mean   :6.494   Mean   : 26.11   Mean   :-74.83  
##  3rd Qu.:71.89   3rd Qu.:6.588   3rd Qu.: 40.02   3rd Qu.:-74.83  
##  Max.   :99.40   Max.   :8.000   Max.   :126.69   Max.   :-74.82  
##       lat                  geometry 
##  Min.   :-8.758   POINT        :24  
##  1st Qu.:-8.758   epsg:4326    : 0  
##  Median :-8.757   +proj=long...: 0  
##  Mean   :-8.757                     
##  3rd Qu.:-8.756                     
##  Max.   :-8.754

Bubbgle Map:

Honoria_elementos4 <- Honoria_elementos3 %>% st_set_geometry(NULL)
coordinates(Honoria_elementos4) <- ~Este+Norte
class(Honoria_elementos4)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

summary(Honoria_elementos4)

## Object of class SpatialPointsDataFrame
## Coordinates:
##           min     max
## Este   518642  519556
## Norte 9031868 9032343
## Is projected: NA 
## proj4string : [NA]
## Number of points: 24
## Data attributes:
##        Cd              Pb              Zn              pH       
##  Min.   :1.100   Min.   :4.740   Min.   :11.80   Min.   :5.410  
##  1st Qu.:1.350   1st Qu.:5.247   1st Qu.:44.79   1st Qu.:6.160  
##  Median :1.575   Median :5.850   Median :58.70   Median :6.375  
##  Mean   :1.634   Mean   :5.929   Mean   :58.66   Mean   :6.494  
##  3rd Qu.:1.863   3rd Qu.:6.525   3rd Qu.:71.89   3rd Qu.:6.588  
##  Max.   :2.550   Max.   :7.240   Max.   :99.40   Max.   :8.000  
##        Ce              lng              lat        
##  Min.   :  4.50   Min.   :-74.83   Min.   :-8.758  
##  1st Qu.:  9.00   1st Qu.:-74.83   1st Qu.:-8.758  
##  Median : 13.50   Median :-74.83   Median :-8.757  
##  Mean   : 26.11   Mean   :-74.83   Mean   :-8.757  
##  3rd Qu.: 40.02   3rd Qu.:-74.83   3rd Qu.:-8.756  
##  Max.   :126.69   Max.   :-74.82   Max.   :-8.754

bubble(Honoria_elementos4, "Cd", col = c("green","green"), main = "Cadmium Concentrations (ppm)")

Grilla generada por Meuse Data:

data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))

# Hacemos una grilla regular con `SpatialPolygonsDataFrame`
proj4string(Honoria_elementos4) <- CRS("+proj=utm +zone=18 +south +datum=WGS84") #definimos proyección

Honoria_elementos4  <- "+proj=longlat +datum=WGS84 +no_defs"

grd <- makegrid(Honoria_elementos4, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos


# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
  coords      = grd, 
  proj4string = CRS(proj4string(Honoria_elementos4))
)


# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts))) +
  geom_point(aes(x, y))

gridded(grd_pts) <- TRUE
grd_pts <- as(grd_pts, "SpatialPixels")

Limite2 <- readOGR(dsn="Data_Suelos/shp/Cacao_shp/Limite_Cacao..shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\CursoDAG\DataAnalisisGeociencias\Modulos\ModuloII\Data_Suelos\shp\Cacao_shp\Limite_Cacao..shp", layer: "Limite_Cacao."
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  Id

grd <- makegrid(Limite2, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos

# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
  coords      = grd, 
  proj4string = CRS(proj4string(Limite2))
)

# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[Limite2, ]

grd_pts_in2 <- grd_pts_in
proj4string(grd_pts_in2) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
gridded(grd_pts_in2) <- TRUE
grd_pts_in2 <- as(grd_pts_in2, "SpatialPixels")

Honoria_elementos5 <- Honoria_elementos4
proj4string(Honoria_elementos5) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
cd.idw = idw(Cd~1, Honoria_elementos5, grd_pts_in)

## [inverse distance weighted interpolation]

spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")

cd.idw = idw(Cd~1, Honoria_elementos5, grd_pts_in2)

## [inverse distance weighted interpolation]

spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")

m <- vgm(.59, "Sph", 874, .04) #setear bien
x <- krige(log(Cd)~1, Honoria_elementos5, grd_pts_in2, model = m)

## [using ordinary kriging]

spplot(x["var1.pred"], main = "ordinary kriging predictions")

spplot(x["var1.var"],  main = "ordinary kriging variance")

Conclusiones

• El analisis de informacion se puede realizar casi completamente en el software R y Rstudio, lo que permite automatizar el proceso en caso de estudios similares.

• El análisis realizado es el resultado de un estudio exitoso en una revista indexada de nivel Q1 en Heavy Metal Pollution and Its Effects on Agriculture en MDPI, El link es el siguiente: Agronomy MDPI.